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ABSTRACT 

Wc present a five-wave Riemann solver for the equations of ideal relativistic magneto- 
hydrodynamics. Our solver can be regarded as a relativistic extension of the five- wave 
HLLD Riemann solver initially developed by Miyoshi and Kusano for the equations 
of ideal MHD. The solution to the Riemann problem is approximated by a five wave 
pattern, comprised of two outermost fast shocks, two rotational discontinuities and a 
contact surface in the middle. The proposed scheme is considerably more elaborate 
than in the classical case since the normal velocity is no longer constant across the ro- 
tational modes. Still, proper closure to the Rankine-Hugoniot jump conditions can be 
attained by solving a nonlinear scalar equation in the total pressure variable which, 
for the chosen configuration, has to be constant over the whole Riemann fan. The 
accuracy of the new Riemann solver is validated against one dimensional tests and 
multidimensional applications. It is shown that our new solver considerably improves 
over the popular HLL solver or the recently proposed HLLC schemes. 
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. 1 MOTIVATIONS 

1 Relativistic flows are involved in many of the high-energy 
' astrophysical phenomena, such as, for example, jets in ex- 
, tragalactic radio sources, accretion flows around compact 
■ objects, pulsar winds and 7 ray bursts. In many instances 
the presence of a magnetic field is also an essential ingredient 
. for explaining the physics of these objects and interpreting 
' their observational appearance. 

Theoretical understanding of relativistic phenomena is 
subdue to the solution of the relativistic magnetohydrody- 
namics (RMHD) equations which, owing to their high degree 
of nonlinearity, can hardly be solved by analytical meth- 
ods. For this reason, the modeling of such phenomena has 
prompted the search for efflcient and accurate num erical 
formu lations. In this respect, Godunov-type schemes (|Tord 
1 19971 ) have gained increasing popularity due to their ability 
and robustness in accurately describing sharp flow disconti- 
nuities such as shocks or tangential waves. 

One of the fundamental ingredient of such schemes is 
the exact or approximate solution to the Riemann problem, 
i.e., the decay between two constant states separated by a 
discont inuity. Unfortun ately the use of an exact Riemann 
solver (|Giacomazzo fc RczzoUa 200 6) is prohibitive because 
of the huge computational cost related to the high degree of 
nonUnearities present in the equations. Instead, approximate 
methods of solution are preferred. 
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Linearized solv ers (|Komissarovl 1 19991 : iBalsaral I2OOII : 
iKoldoba et "aLll2002l ) rely on the rather convoluted eigenvec- 
tor decomposition of the underlying equations and may be 
prone to numerical pathologies leading to negative density 
or pressures inside the solution (|Einfeldt et al.|[l991^ . 

Characteristic-free algorithms based on the Ru- 
sanov Lax-Friedri chs or the Harten-Lax-van Leer (HLL, 
iHarten et "allll983l ) formulations are sometime preferred due 
to their ease of implementation and positivity properties. 
Implementation of such algorithm s can be found in the 
codes d escribed by iGammie et al.l (l2003h: iLeism ann ct al.l 
(j2005 ): iDerZanna et all (|2007l ): Ivan der Holst'eTal.. l(2Qim . 
Although simpler, the HLL scheme approximates only two 
out of the seven waves by collapsing the full structure of 
the Riemann fan into a single average state. These solvers, 
therefore, are not able to resolve intermediate waves such as 
Alfven, contact and slow discontinuities. 

Attempts to restore the middle contact (or en- 
tropy) wave (HLLC, i n itially devised for the Euler 
equations by pToro et al] Il994h have been proposed by 
iMignone et al.l (|2005|) in the case o f purely transversal 
flelds and by iMignone fc Bodd (|2006l ) (MB from now on), 
iHonkkila fc JanhunenI (|2007^ in the more general case. 
These schem es prov i de a r elativistic exten sion of the work 
proposed by ICurskil l|2004 l and Q (|2005h for the classical 
MHD equations. 

HLLC solvers for the equations of MHD and RMHD, 
however, still do not capture slow discontinuities and Alfven 
waves. Besides, direct application of the HLLC solver of MB 
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to genuinely 3D problems was shown to suffer from a (poten- 
tial) pathological singularity when the component of mag- 
netic field normal to a zone interface approaches zero. 

A step forward in reso lving intermediate wave struc- 
tures was then performed bv lMivoshi &i Kusand (|2005l ') (MK 
from now on) who, in the context of Newtonian MHD, in- 
troduced a four state solver (HLLD) restoring the rotational 
(Alfven) discontinuities. In this paper we propose a gener- 
alization of Miyoshi & Kusano approach to the equations 
of relativistic MHD. As we shall see, this task is greatly 
entangled by the different nature of relativistic rotational 
waves across which the velocity component normal to the 
interface is no longer constant. The proposed algorithm has 
been implemen ted in the PLUTO code for astrophysical 
fluid dynamics (jMignone et al.ll2007l ) which embeds a com- 
putational infrastructure for the solution of different sets of 
equations (e.g., Euler, MHD or relativistic MHD conserva- 
tion laws) in the finite volume formalism. 

The paper is structured as follows: in ij2]we briefly re- 
view the equations of relativistic MHD (RMHD) and formu- 
late the problem. In ^the new Riemann solver is derived. 
Numerical tests and astrophysical applications are presented 
in 34] and conclusions are drawn in SJ5] 



2 BASIC EQUATIONS 

The equations of relativistic magnetohydrodynamics 
(RMHD) are derived under the physical assumptions of 
constant magnetic permeability and infinite conductivity , 
appropriate for a p erfectly conducting fluid (jAnild 1 19891 : 
iLichnerowicd Il967h . In divergence form, they express 
particle number and energy-momentum conservation: 



9m i.P^'' 



[wg + b^) u^u - h^V" + ( Pg + — ) 7? 



0, (2) 



d^iu^b" -^u"b'') = 0, (3) 

where p is the rest mass density, u'^ = 7(1,1;) is the four- 
velocity (7 = Lorentz factor, v = three velocity), lOg and Pg 
are the gas enthalpy and thermal pressure, respectively. The 
covariant magnetic field b'^ is orthogonal to the fiuid four- 
velocity {u''bfj_ = 0) and is related to the local rest frame 
field B by 



71; ■ B, \- J {v ■ B) V 

7 



(4) 

B)^ is the squared 



In Eq. b^ = b^'b^, = BVt^ + (« ' 
magnitude of the magnetic field. 

The set of equations (Il|-(|3]) must be complemented by 
an equation of state which may be taken as the constant 
r-law: 

r 

rPa > (5) 



P 



V 



where F is the specific heat ratio. Alternative equa tions of 
state (see, for example. iMignone fc McKinnevl[2007l 'l may be 
adopted. 

In the following we will be dealing with the one dimen- 
sional conservation law 



dU ^ 9F _^ 

dt dx ' 



(6) 



which follows directly from Eq. (IT])-(|3]) by discarding contri- 
butions from y and 2. Conserved variables and correspond- 
ing fiuxes take the form: 



U 



( ^ \ 

E 

\ B" J 



\ 



(7) 



where k = x,y, z, D = p"/ is the the density as seen from 
the observer's frame while, introducing w = Wg + b^ (total 
enthalpy) and p = pg + &^/2 (total pressure), 
uOik TP 0.0 ,0,0 



E = wu u —00 — : 



(8) 



are the momentum and energy densities, respectively. 5kx is 
the Kronecker delta symbol. 

Note that, since Fsx = 0, the normal component of 
magnetic field (B^) does not change during the evolution 
and can be regarded as a parameter. This is a direct conse- 
quence of the V • B = condition. 

A conservative discretization of Eq. ((6)1 over a time step 
At yields 



At 



(9) 



Ax V '+5 •' »-f 

where Ax is the mesh spacing and f^^i is the upwind nu- 
merical flux computed at zone faces x^_^_i by solving, for 
t" < t < the initial value problem defined by Eq. (O 

together with the initial condition 



0, (1) u{x,n^ 



Ui 



for X < x^_ 
for X > X- 



(10) 



where Ul and Ur are discontinuous left and right constant 
states on either side of the interface. This is also known as 
the Riemann problem. For a first order scheme, Ul ~ Ui 
and Ur = Ui+i. 

The decay of the initial discontinuity given by Eq. (|10|l 
leads to the formation of a self-similar wave pattern in the 
X — t plane where fast, slow, Alfven and contact modes can 
develop. At the double end of the Riemann fan, two fast 
magneto-sonic waves bound the emerging pattern enclosing 
two rotational (Alfven) discontinuities, two slow magneto- 
sonic waves and a contact surface in the middle. The same 
patterns is also found in classical MHD. Fast and slow 
magneto-sonic disturbances can be either shocks or rar- 
efaction waves, depending on the pressure jump and the 
norm of the magnetic field. All variables (i.e. density, ve- 
locity, magnetic field and pressure) change discontinuously 
across a fast or a slow shock, whereas thermodynamic quan- 
tities such as thermal pressure and rest density remain con- 
tinuous when crossing a relativistic Alfven wave. Contrary 
to its classical counterpart, however, the tangential compo- 
nents of magnetic field trace ellipses instead of circles and 
the normal component of the velocity is no longer contin- 
uous across a rotational discontinuity, iKomissaro 3 (|l99'if ). 
Finally, through the contact mode, only density exhibits a 
jump while thermal pressure, velocity and magnetic field re- 
main continuous. 
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Figure 1. Approximate structure of the Riemann fan introduced 
by the HLLD solver. Tlie initial states and Un arc connected 
to each other through a set of five waves representing, clockwise, a 
fast shock A^, a rotational discontinuity A^Li a- contact wave Ac, a 
rotational discontinuity A^^ and a fast shock Aij. The outermost 
states, and Un are given as input to the problem, whereas 
the others must be determined consistently solving the Rankine- 
Hugoniot jump conditions. 



The complete analytical solution to the Riemann prob- 
le m in RMHD has been r e cently derived in closed form 
by iGiacomazzo fc Rezzollal (|2006h and number of proper- 
ties regardinp; s i mple waves are als o well established, see 
lAnile fc Pennisil (| 19871 ): I Anil j (|l989l ). 

For the special case in which the component of the mag- 
netic field normal to a zone interface vanishes, a degeneracy 
occurs where tangential, Alfven and slow waves all propa- 
gate at the speed of the fluid and the s o lution simplifies to 
a three- wave pattern, see iRomero et al.1 (|2005l ). 

The high degree of nonlinearity inherent to the RMHD 
equations makes seeking for an exact solution prohibitive in 
terms of computational costs and efficiency. For this reasons, 
approximate methods of solution are preferred instead. 



3 THE HLLD APPROXIMATE RIEMANN 
SOLVER 

Without loss of generality, we place the initial discontinuity 
at a; = and set t" = 0. 

Following MK, we make the assumption that the Rie- 
mann fan can be divided by 5 waves: two outermost fast 
shocks, Xr and Xl, enclosing two rotational discontinuities, 
XaL and XaR, separated by the entropy (or contact) mode 
with speed Ac. Note that slow modes are not considered in 
the solution. The five waves divide the x — t plane into the 
six regions shown in Fig[T] corresponding (from left to right) 
to the 6 states Ua with a — L, aL, cL, cR, aR, R. 

The outermost states {U l and U r) are given as input 
to the problem, while the remaining ones have to be deter- 
mined. In the typical approach used to construct HLL-based 
solvers, the outermost velocities Xl and Xr are also provided 
as estimates from the input left and right st ates. As in M B, 
we choose to use the simple Davis estimate (|Davisll 19881 ). 

Across any given wave A, states and fluxes must satisfy 
the jump conditions 



where -f and — identify, respectively, the state immedi- 
ately ahead or behind the wave front. Note that for consis- 
tency with the integral form of the conservation law over 
the rectangle [A_l At, AijAt] x [0, At] one has, in general. 
Fa 7^ FiUci), except of course for a — L or a = R. 

Across the fast waves, we will make frequent use of 



Rl ~ XlU L — F I 



Rr = XrUr — Fr . 



(12) 



V 
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= 0, 


P 


= 




Ac 











which are known vectors readily obtained from the left and 
right input states. A particular component of R is selected 
by mean of a subscript, e.g., Rd is the density component 
oiR. 

A consistent solution to the problem has to satisfy the 
7 nonlinear relations implied by Eq. (|ll|l for each of the 5 
waves considered, thus giving a total of 35 equations. More- 
over, physically relevant solutions must fulfill a number of 
requirements in order to refiect the characteristic nature of 
the considered waves. For this reason, across the contact 
mode, we demand that velocity, magnetic field and total 
pressure be continuous: 



(13) 



and require that Ac = «c , i-e-, that the contact wave moves 
at the speed of the fluid. However, density, energy and total 
enthalpy may be discontiimous. On the other hand, through 
the rotational waves XaL and XaR, scalar quantities such as 
total pressure and enthalpy are invariant whereas all vector 
components (except for B^) experience jumps. 

Since slow magnetosonic waves are not considered, we 
naturally conclude that only the total pressure remains con- 
stant throughout the fan, contrary to Newtonian MHD, 
where also the velocity normal to the interface {v^) is left 
unchanged across the waves. This is an obvious consequence 
of the different nature of relativistic Alfven waves across 
which vector fields like u** and trace ellipses rather than 
circles. As a consequence, the normal component of the ve- 
locity, , is no longer invariant in RMHD but experiences 
a jump. These considerations along with the higher level of 
complexity of the relativistic equations makes the extension 
of the multi-state HLL solver to RMHD considerably more 
elaborate. 

Our strategy of solution is briefly summarized. For each 
state we introduce a set of 8 independent unknowns: V — 
{D,v'^ , B^ , B" jWjp} and write conservative variables 
and fluxes given by Eq. (0 as 



Ua 



D 

wy^v'' - b°b^ 



\ 



\ 



(14) 



-^p-b"b" 
( \ 

where k = x,y,z labels the vector component, a is the state 



Fa 



\ 



(15) 
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and is computed directly from Q. We proceed by solv- 
ing, as function of the total pressure p, the jump conditions 

across the outermost waves Xl and \r. By requiring 
that total pressure and Alfven velocity do not change across 
each rotational modes, we find a set of invariant quantities 
across XaL and XaR- Using these invariants, we express states 
and fluxes on either side of the contact mode (a = cL, cR) in 
terms of the total pressure unknown only. Imposing continu- 
ity of normal velocity, v^i^{p) = v^ii{p), leads to a nonlinear 
scalar equation in p, whose zero gives the desired solution. 

Once p has been found to some relative accuracy (typi- 
cally 10~^), the full solution to the problem can be written 
as 



Fl 




if 


Xl>0 


FaL 




if 


Al < < XaL 


FaL + Ac {U cL 


-UaL) 


if 


XaL <0 < Ac 


FaR + Ac {U cR 


-UaR) 


if 


Ac < < XaR 


FaR 




if 


XaR <0 < Xr 


Fr 




if 


Xr<0 



(16) 



where UaL, UaR. are computed in ii3.ll UcL,UcR in ^3.31 
and Fa — F + Xa{Ua — U) (for a — aL or a = aR) follow 
from the jump conditions. The wave speeds XaL, XaR and Ac 
are computed during the solution process. 

Here and in what follows we adopt the convention that 
single subscripts like a (or c) refers indifferently to aL, aR (or 
cL, cR). Thus an expression like Wc = Wa means Wcl = WaL 
and WcR — WaR- 



3.1 Jump Conditions Across the Fast Waves 

We start by explicitly writing the jump conditions across 
the outermost fast waves: 

(A-i;")L> = Rd, (17) 
[X-v^w-f^" +b'' {b^ -Xb") -p5,^ = , (18) 

(A-w")«;7'-Ap + fe"(6"-A6'') = Re, 
(A-«")B'' 



(19) 
(20) 



where, to avoid cluttered notations, we omit in this section 
the Q = aL (when A = AL)ora = aR (when A = Xr) index 
from the quantities appearing on the left hand side. Simi- 
larly, the R's appearing on the right hand sides of equations 
H17p - (|20|) are understood as the components of the vector 
Rl (when A = Al) or Rr (when A = Xr), defined by Eq. 

The jump conditions of Faraday's law allow to express 
the magnetic field as a function of velocities alone, 

k = y,z. (21) 



X — v^' 



for 



The energy and momentum equations can be combined to- 
gether to provide an explicit functional relation between the 
three components of velocity and the total pressure p. To 
this purpose, we first multiply the energy equation (|19p 
times v'' and then subtract the resulting expression from 
the jump condition for the fc-th component of momentum. 



Eq. IISI). Using Eq. to get rid of the term, one finds 
after some algebra: 

B" (b" -Rb-v^ - - Xv"^ = - v^Re , (22) 

with defined by (|2H) . The system can be solved for v'^ 
giving 

, _ B^AB'' + XC)- (A + G){p + Rrn^) 



X 

y ^ QRmv + Rbv [C + B" {XR^. - Re)] 

2 _ QRm' + Rb' [C + B^ (^XRm" — Re)] 

" ^ X 
where 

A = ~XRe+p{1-X^) , 

G = Rbv Rbv -\- Rb^Rb^ 1 

C = RmV Rbv + Rju^Rb^ 1 

Q = -A-G+(B")' (l- A^) 

X = B"" {AXB'' + C) - {A + G) {Xp + Re) 



(23) 
(24) 
(25) 

(26) 
(27) 
(28) 
(29) 
(30) 



Once the velocity components are expressed as func- 
tions of p, the magnetic field is readily found from (|2H) . 
while the total enthalpy can be found using its definition, 
w = {E + p)/'y'^ + {v ■ B)^, or by subtracting Re from the 
inner product v'' ■ Rm, giving 

Re — V ■ Rm 



w — p + 



X — 



(31) 



where Rm = {Rm^ , Rmv , Rtn')- Although equivalent, we 
choose to use this second expression. Since the are func- 
tions of p alone, the total enthalpy w is also a function of 
the total pressure. 

The remaining conserved quantities in the a = aL or 
a = aR regions can be computed once p has been found: 

D = -r^, (32) 



E = 



X — 

Re + pv"' -jv-B) B" 
X — 



m'= = {E + p)v'' -{v-B)B'' . 



(33) 
(34) 



One can verify by direct substitution that the previous equa- 
tions together with the corresponding fiuxes, Eq. (|15|) . sat- 
isfy the jump conditions given by (|17p - (|20|l . 



3.2 Jump Conditions across the Alfven waves 

Across the rotational waves one could, in principle, proceed 
as for the outer waves, i.e., by explicitly writing the jump 
conditions. However, as we shall see, the treatment greatly 
simplifies if one introduces the four vector 



(7^ = Tju'^ + b^ 



with T) — ±sign(_B'^)\A() 



(35) 



where, for reasons that will be clear later, we take the plus 
(minus) sign for the right (left) state. From a'^ we define the 
spatial vector K = {K'^ , , K") with components given by 



_ — * + 

- ^ + 7aO 



(36) 
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The vector K has some attractive properties, the most re- 
markable of which is that the x compone nt coincides with 
the propagation speed of the Alfven wave (|Anilelll989l ). For 
this reason, we are motivated to define Aa = K^, where 
the subscript a stands for either the left or right rota- 
tional wave (i.e. aL or aR) since we require that both K'^ 
and p are invariant across the rotational discontinuity, i.e., 

— = pc — Pa ~ 0, a, property certainly shared by the 
exact solution. As we will show, this choice naturally reduces 
to the expressions found by MK in the non-relativistic limit. 

Indeed, setting Aa = Kl — and using Eq. (|36} to 
express as functions of K'^ , the jump conditions simplify 
to 



DB^ 



- p5kjL 







= 



= 0, 



(37) 
(38) 
(39) 
(40) 



Since also [p]a„ ~ 0, the previous equations further imply 
that (when / 0) also D/ija^), w, and K"" do not 
change across Aa: 



KaL — KcL '- 
KaR = KcR. 



Kl, 

■Kr, 



TjaL = rich — TjL 
rjaR = VcR = VR 



(41) 
(42) 



Being invariant, K can be computed from the state 
lying to the left (for Aai) or to the right (for XaR) of the 
discontinuity, thus being a function of the total pressure p 
alone. Instead of using Eq. (|36|l . an alternative and more 
convenient expression may be found by properly replacing 
n*' with K'' in Eq. (fT7)l - (PIJ)) . After some algebra one finds 
the simpler expression 

jy^k ^ Rmk +pSkx + RskV . 

Xp + Re + B-t^ ' ^^^> 

still being a function of the total pressure p. 

Note that, similarly to its non relativistic limit, we can- 
not use the equations in (|37p - (|40[) to compute the solu- 
tion across the rotational waves, since they do not provide 
enough independent relations. Instead, a solution may be 
found by considering the jump conditions across both rota- 
tional discontinuities and properly matching them using the 
conditions at the contact mode. 



one has = = B^, where 

[b'=(A - u") + B^'v''] - [b'=(A - L.") + B'u'-'] 



B. 



XaR — XaL 

Since quantities in the aL and aR regions are given in terms 
of the p unknown, Eq. (|45|l are also functions of p alone. 

At this point, we take advantage of the fact that cr'^it^ = 
—77 to replace ja" with 77/(1 — AT • w) and then rewrite (|36|l 
as 



(45) 



B' 



{l-K -v) 



for 



x,y,z. 



(46) 



The previous equations form a linear system in the velocity 
components and can be easily inverted to the left and to 
the right of the CD to yield 



B^jl-K') 
r]~K B 



for 



x,y,z. 



(47) 



which also depend on the total pressure variable only, with 
w and K'^ being given by (|3T]| and (|43)) and the B^ 's being 
computed from Eq. (|45|) . Imposing continuity of the normal 
velocity across the CD, v^j^ — = 0, results in 



where 



Ys{p) 



Yl 



0, 



77sAi4'- -Ks 



S — L, R . 



(48) 



(49) 



is a function of p only and Be = AK^Bc is the numerator 
of (fiS)) and AJf^ = K^j^ — K^j^. Equation (gSI is a nonlinear 
function in p and must be solved numerically. 

Once the iteration process has been completed and p 
has been found to some level of accuracy, the remaining 
conserved variables to the left and to the right of the CD 
are computed from the jump conditions across XaL and XaR 
and the definition of the flux, Eq. p5|) . Specifically one has, 
for {c = cL, a = aL} or {c = cR, a = aR}, 



Da 



Xa - Vl 



Aa - 

XgEa -ml+ pv^ ~ {vc ■ Be) B^ 

Xa —Vc 

ml = {Ec + p)v'l ~ {vc ■ Be) B^ . 

This concludes the derivation of our Riemann solver. 



(50) 

(51) 
(52) 



3.3 Jump Conditions across the Contact wave 

At the contact discontinuity (CD) only density and total 
enthalpy can be discontinuous, while total pressure, normal 
and tangential fields are continuous as expressed by Eq. p3p . 

Since the magnetic field is a conserved quantity, one 
can immediately use the consistency condition between the 
innermost waves AaL and XaR to find B'' across the CD. 
Indeed, from 

iXe^XaL)UcL+{XaR-Xe)UcR= (44) 
= XaRUaR — XahU aL — F aR + F aL 



3.4 Full Solution 

In the previous sections we have shown that the whole set 
of jump conditions can be brought down to the solution of 
a single nonlinear equation, given by (|48|l . in the total pres- 
sure variable p. In the particular case of vanishing normal 
component of the magnetic field, i.e. B^ —> 0, this equation 
can be solved exactly as discussed in i]3.4.1l 

For the more general case, the solution has to be found 
numerically using an iterative method where, starting from 
an initial guess , each iteration consists of the following 
steps: 
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• given a new guess value p''"' to the total pressure, start 
from Eq. (|23|) - (|25|l to express VaL and VaR as functions of 
the total pressure. Also, express magnetic fields BaL, BaR 
and total enthalpies wl, wr using Eq. pif) and Eq. (|3H) . 
respectively. 

• Compute KaL and KaR using Eq. (|43|l and the trans- 
verse components of Be using Eq. (|45|l . 

• Use Eq. (|48|) to find the next improved iteration value. 

For the sake of assessing the validity of our new solver, we 
choose the secant method as our root-finding algorithm. The 
initial guess is provided using the following prescription: 



when (B^)Vp''" < 0.1, 
otherwise , 



(53) 



where pi^" is the total pressure computed from the HLL av- 
erage state whereas po is the solution in the — limiting 
case. Extensive numerical testing has shown that the total 
pressure p*^" computed from the HLL average state provides, 
in most cases, a sufficiently close guess to the correct physi- 
cal solution, so that no more than 5 — 6 iterations (for zones 
with steep gradients) were required to achieve a relative ac- 
curacy of 10~^. 

The computational cost depends on the simulation set- 
ting since the average number of iterations can vary from 
one problem to another. However, based on the results pre- 
sented in ^ we have found that HLLD was at most a factor 
of ~ 2 slower than HLL. 

For a solution to be physically consistent and well- 
behaved, we demand that 



WR> p, vIj^KXr, V^h < XaR : 



(54) 



hold simultaneously. These conditions guarantee positivity 
of density and that the correct eigenvalue ordering is always 
respected. We warn the reader that equation (|48p may have, 
in general, more than one solution and that the conditions 
given by (|54|l may actually prove helpful in selecting the cor- 
rect one. However, the intrinsic nonlinear complexity of the 
RMHD equations makes rather arduous and challenging to 
prove, a priori, both the existence and the uniqueness of a 
physically relevant solution, in the sense provided by (|54|) . 
On the contrary, we encountered sporadic situations where 
none of the zeroes of Eq. (|48|l is physically admissible. For- 
tunately, these situations turn out to be rare eventualities 
caused either by a large jump between left and right states 
(as at the beginning of integration) or by under- or over- esti- 
mating the propagation speeds of the outermost fast waves, 
Xl and Xr. The latter conclusion is supported by the fact 
that, enlarging one or both wave speeds, led to a perfectly 
smooth and unique solution. 

Therefore, we propose a safety mechanism whereby we 
switch to the simpler HLL Riemann solver whenever at least 
one or more of the conditions in (|54|l is not fulfilled. From 
several numerical tests, including the ones shown here, we 
found the occurrence of these anomalies to be limited to few 
zones of the computational domain, usually less than 0.1% 
in the tests presented here. 

We conclude this section by noting that other more 
sophisticated algorithms may in principle be sought. One 
could, for instance, provide a better guess to the outer wave- 



speeds Xl and Xr or even modify them accordingly until a 
solution is guaranteed to exist. Another, perhaps more use- 
ful, possibility is to bracket the solution inside a closed in- 
terval [pmin,Pmax] where Pmin and Pmax may be found from 
the conditi ons (1541). Using a n alternative root finder, such 
as Kidder (|Press et al.lll993 ). guarantees that the solution 
never jump outside the interval. However, due to the small 
number of failures usually encountered, we do not think 
these alternatives could lead to a significant gain in accu- 
racy. 

3.4-1 Zero normal field limit 

In the limit B^ — + a degeneracy occurs where the Alfven 
(and slow) waves propagate at the speed of the contact mode 
which thus becomes a tangential discontinuity. Across this 
degenerate front, only normal velocity and total pressure re- 
main continuous, whereas tangential vector fields are subject 
to jumps. 

This case does not pose any serious difficulty in our 
derivation and can be solved exactly. Indeed, by setting 
B"" = in Eq. (|43)) and (|48)) . one immediately finds that 

= = Wc leading to the following quadratic equation 
for p: 



p^ + (ij"' - 



J. „ 



\)p + 



m tp- 



0, 



(55) 



where the superscript "hll" refers to the HLL average state 
or flux given by Eq. (28) or (31) of MB. We note that equa- 
tion j55|_com£idBswit^ the derivation given by MB (see 
also iMignone et~al]|2005^ in the same degenerate case and 
the positive root gives the correct physical solution. The in- 
termediate states, U cL and U cR, loose their physical mean- 
ing as B^ but they never enter the solution since, as 
XaL, ^aR — » Ac, Only U aL U aR. will havc a uouzcro finite 
width, see Fig. [T] 

Given the initial guess, Eq. (|53p . our proposed approach 
does not have to deal separately with the B^ ^ a nd B^ — 
cases (as in MB and lHonkkila fc Janhunenll2007l ) and thus 
solves the issue raised by MB. 

3.4-2 Newtonian Limit 

We now show that our derivation reduces to the HLLD 
Riemann solver found by MK under the appropriate non- 
relativistic limit. We begin by noticing that, for v/c^ 0, the 
velocity and induction four-vectors reduce to (li'^'^) 
and 6^ (0,5'°), respectively. Also, note that Wg,w p 
in the non-relativistic limit so that 

B^ 



K 



V + s 



VP 



(56) 



and thus cannot change across Xa- Replacing (|17|l - (|18[ l 
with their non-relativistic expressions and demanding — 
gives, in our notations, the following expressions: 



P = (B 



Rr,d — Rl,d ' 

RL,m^ RRji — RR,m^ Rr^ 



Rf 



Rl.d 



(57) 
(58) 



which can be shown to be identical to Eqns (38) and (41) 
of MK. With little algebra, one can also show that the re- 
maining variables in the aL and aR regions reduce to the 
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corresponding non-relativistic expressions of MK. Similarly, 
the jump across the rotational waves are solved exactly in 
the same fashion, that is, by solving the integral conservation 
laws over the Riemann fan. For instance, Eq. (|45p reduces 
to equation (61) and (62) of MK. These results should not 
be surprising since, our set of parameters to write conserved 
variables and fluxes is identical to the one used by MK. The 
only exception is the energy, which is actually written in 
terms of the total enthalpy. 



4 NUMERICAL TESTS 

We now evaluate, in ij4.1l the accuracy of the proposed 
HLLD Riemann solver by means of selected one dimen- 
sional shock tube problems. Applications of the solver to 
multi-dimensional problems of astrophysical relevance are 
presented in i]4.2l 



4.1 One Dimensional Shock Tubes 

The initial condition is given by Eq. (|10|) with left and right 
states defined by the primitive variables listed in Table 1. 
The computational domain is chosen to be the interval [0, 1] 
and the discontinuity is placed at a; = 0.5. The resolution 
Nx and final integration time can be found in the last two 
columns of Table 1. Unless otherwise stated, we employ the 
constant F— law with F — 5/3. The RMHD equations are 
solved using the first-order accurate scheme ([UJ with a CFL 
number of 0.8. 

Numerical results are compared to the HLLC Riemann 
solver of MB and to the simpler HLL scheme and the accu- 
racy is quantified by computing discrete errors in L-1 norm: 



ELI = 



— Qi Axi , 



(59) 



where qi is the first-order numerical solution (density or 
magnetic field), ql"^ is the reference solution at Xi and Axi 
is the mesh spacing. For tests 1,2,4 we obtained a refer- 
ence solution using the second-order scheme of MB on 3200 
zones and adaptive mesh refinement with 6 levels of re- 
finement (equivalent resolution 204,800 grid points). Grid 
adaptivity in one dimension has been incorporated in the 
PLUT O code using a block-s tructured grid approach fol- 
lowing HergeLTColini (|l989l ). For test 3, we use the exact 
numerical solution available from iGiacomazzo fc Rezzollal 
(|2006l ). Errors (in percent) are shown in Fig. 1111 



4-1.1 Exact Resolution of Contact and Alfven 
Discontinuities 

We now show that our HLLD solver can capture exactly 
isolated contact and rotational discontinuities. The initial 
conditions are listed at the beginning of Table 1. 

In the case of an isolated stationary contact wave, only 
density is discontinuous across the interface. The left panel 
in Fig. [2] shows the results at t = 1 computed with the 
HLLD, HLLC and HLL solvers: as expected our HLLD pro- 
duces no smearing of the discontinuity (as does HLLC). On 
the contrary, the initial jump broadens over several grid zone 
when employing the HLL scheme. 



Isolated 
Contact 
Wove 



• HLLD 
XHLLC 
+ HLL 



0.6 0.8 



Isolated < 

Rotational 

Wave 



• HLLD 
XHLLC 
+ HLL 



0.2 L_ 
0.0 



0.5 0.8 1.0 



Figure 2. Results for the isolated contact (left panel) and rota- 
tional (right panel) waves at t = 1. Density and y component of 
magnetic field are shown, respectively. The difl'erent symbols show 
results computed with the new HLLD solver (filled circles), the 
HLLC solver (crosses) and the simpler HLL solver (plus signs). 
Note that only HLLD is able to capture exactly both discontinu- 
ities by keeping them perfectly sharp without producing any grid 
difi^usion effect. HLLC can capture the contact wave but not the 
rotational discontinuity, whereas HLL spreads both of them on 
several grid zones. 







1.0 0.0 0-2 0.4 0.6 0.8 1.0 0.0 0.2 0,4 0.6 0.8 1,0 



Figure 3. Relativistic Brio-Wu shock tube test at t = 0.4. Com- 
putations are carried on 400 zones using the HLLD (solid line), 
HLLC (dashed line) and HLL (dotted line) Riemann solver, re- 
spectively. The top panel shows, from left to right, the rest mass 
density, gas pressure, total pressure. The bottom panel shows the 
X and y components of velocity and the y component of magnetic 
field. 



Across a rotational discontinuity, scalar quantities such 
as proper density, pressure and total enthalpy are invariant 
but vector fields experience jumps. The left and right states 
on either side of an exact rotational discontinuity can be 
found using the procedure outlined in the Appendix. The 
right panel in Fig. [2] shows that only HLLD can successfully 
keep a sharp resolution of the discontinuity, whereas both 
HLLC and HLL spread the jump over several grid points 
because of the larger numerical viscosity. 



4.1.2 Shock Tube 1 

The first shock tube test is a relativistic exten sion of the Brio 
Wu magnetic shock tube (|Brio fc "Wulll988h and has also 
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Table 1. Initial conditions for the test problems discussed in the text. The last two columns give, respectively, the final integration time 
and the number of computational zones used in the computation. 










0.0 






— ?r 


2.5 






k 










2-0 






0.4 


\ 






1.5 




















□.6 






















1.0 



0.50 0.55 0.60 0.65 0.70 



0.50 0.55 0.6O 0.65 0.70 



0.50 0.55 0.60 0.65 0.70 



Figure 4. Enlargement of the central region of Fig. |3] Density 
and the two components of velocity are shown in the left, cen- 
tral and right panels, respectively. Diamonds, crosses and plus 
signs are used for the HLLD, HLLC and HLL Riemann solver, 
respectively. 



been considered bv lBalsaral (j200lf ): lOel Zanna et al.1 (|2003l l 
and in MB. The specific heat ratio is F — 2. The initial 
discontinuity breaks into a left-going fast rarefaction wave, a 
left-going compound wave, a contact discontinuity, a right- 
going slow shock and a right-going fast rarefaction wave. 
Rotational discontinuities are not generated. 

In Figs. |3]|4]we plot the results obtained with the first- 
order scheme and compare them with the HLLC Riemann 
solver of MB and the HLL scheme. Although the resolu- 
tion across the continuous right-going rarefaction wave is 
essentially the same, the HLLD solver offers a considerable 
improvement in accuracy in the structures located in the 
central region of the plots. Indeed, Fig. |4] shows an enlarge- 
ment of the central part of the domain, where the com- 
pound wave (at x ~ 0.51), contact {x ~ 0.6) and slow shock 
{x « 0.68) are clearly visible. Besides the steeper profiles of 
the contact and slow modes, it is interesting to notice that 
the compound wave, composed of a slow shock adjacent to 
a slow rarefaction wave, is noticeably better resolved with 
the HLLD scheme than with the other two. 

These results are supported by the convergence study 
shown in the top left panel of Fig. 1111 demonstrating that the 
errors obtained with our new scheme are smaller than those 
obtained with the HLLC and HLL solvers (respectively). At 
the largest resolution employed, for example, the L-1 norm 
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0.14 . 
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Figure 5. Results for the second shock tube problem at t = 0.55 
on 800 grid points. From left to right, the top panel shows den- 
sity, gas and total pressure. The middle panel shows the three 
components of velocity, whereas in the bottom panel we plot the 
Lorentz factor and the transverse components of magnetic field. 
Solid, dashed and dotted lines are used to identify results com- 
puted with HLLD, HLLC and HLL, respectively. 



errors become ~ 63% and ~ 49% smaller than the HLL and 
HLLC schemes, respectively. 

The CPU times required by the different Riemann 
solvers on this particular test were found to be scale as 
thu '■ thiic : thiid = 1 : 1-2 : 1.9. 



4.1.3 Shock Tube 2 

This test has also been considered in iBalsaral (|200lD and 
in MB and the initial condition comes out as a non-planar 
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Figure 6. Left panel: enlargement of the central region of Fig. 
[5] around the contact wave. Middle and right panels: close-ups 
of the 2 component of velocity and y component of magnetic 
field around the right-going slow shock and Alfven discontinuity. 
Different symbols refer to different Riemann solver, see the legend 
in the left panel. 
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Riemann problem implying that the change in orientation 
of the transverse magnetic field across the discontinuity is 
« 0.557r (thus different from zero or tt). 

The emerging wave pattern consists of a contact wave 
(at X ~ 0.475) separating a left-going fast shock {x « 0.13), 
Alfven wave (x ~ 0.185) and slow rarefaction {x ~ 0.19) 
from a slow shock {x ~ 0.7), Alfven wave {x ~ 0.725) and 
fast shock (x « 0.88) heading to the right. 

Computations carried out with the f^' order accurate 
scheme are shown in Fig. [5] using the HLLD (solid line), 
HLLC (dashed line) and HLL (dotted line). The resolution 
across the outermost fast shocks is essentially the same for 
all Riemann solvers. Across the entropy mode both HLLD 
and HLLC attain a sharper representation of the disconti- 
nuity albeit unphysical undershoots are visible immediately 
ahead of the contact mode. This is best noticed in the the 
left panel of Fig. |B] where an enlargement of the same region 
is displayed. 

On the right side of the domain, the slow shock and 
the rotational wave propagate quite close to each other and 
the first-order scheme can barely distinguish them at a res- 
olution of 800 zones. However, a close-up of the two waves 
(middle and right panel in Fig. [S| shows that the proposed 
scheme is still more accurate than HLLC in resolving both 
fronts. 

On the left hand side, the separation between the Alfven 
and slow rarefaction waves turns out to be even smaller and 
the two modes blur into a single wave because of the large 
numerical viscosity. This result is not surprising since these 
features are, in fact, challenging even for a second-order 
scheme (|Balsarall206ll ). 

Discrete L-1 errors computed using Eq. (|59|l are plot- 
ted as function of the resolution in the top right panel of 
Fig. 111! For this particular test, HLLD and HLLC produce 
comparable errors (~ 1.22% and ~ 1.33% at the highest 
resolution) while HLL performs worse on contact, slow and 
Alfven waves resulting in larger deviations from the refer- 
ence solution. 

The computational costs on 800 grid zones has found 
to be thu ■■ thiic : ihiw = 1 : 1.1 : 1.6. 

4.1.4 Shock Tube 3 

In this test problem we consider the interact ion of two oppo - 
sitely colliding relativis tic streams, see also iBalsaral (|200lD : 
iDel Zanna et al.l (120031 ) and MB. 



Figure 7. Relativistic collision of two oppositely moving streams 
at t = 0.4. From top to bottom, left to right, the panels show 
density p, gas pressure pg, total pressure p, x and y components 
of velocity (i;'^ ) and y component of magnetic field . The 2 
components have been omitted since they are identical to the y 
components. Solid, dashed and dotted lines refer to computations 
obtained with the HLLD, HLLC and HLL solvers. 400 computa- 
tional zones were used in the computations. 







By 





















a.*o 0.45 0.50 0.55 
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Figure 8. Enlargement of the central region in Fig. [7] Filled 
circles crosses and plus sign have the same meaning as in Fig. |6] 
Note the wall heating problem evident in the density profile (left 
panel). Central and right panels show the transverse field profiles. 
Clearly the resolution of the slow shocks (x ^ 0.5it0.07) improves 
from HLL to HLLC and more from HLLC to HLLD. 



After the initial impact, two strong relativistic fast 
shocks propagate outwards symmetrically in opposite direc- 
tion about the impact point, x = 0.5, see Fig. [7] Being a 
co-planar problem (i.e. the initial twist angle between mag- 
netic fields is tt) no rotational mode can actually appear. 
Two slow shocks delimiting a high pressure constant den- 
sity region in the center follow behind. 

Although no contact wave forms, the resolution across 
the slow shocks noticeably improves changing from HLL to 
HLLC and from HLLC to HLLD, see Fig. [7] or the enlarge- 
ment of the central region shown in Fig. |8] The resolution 
across the outermost fast shocks is essentially the same for 
all solvers. 

The spurious density undershoot at the center of the 
grid is a notorious numerical pathology, known as the 
wall heat ing problem, often encountered in Godunov-type 
schemes 



t mg problem, often encountered m (jodunov-type 
(|Nohlll987l : [Ge"hmevr et al.lll997l ). It consists of an 
undesired entropy buildup in a few zones around the point 
of symmetry. Our scheme is obviously no exception as it can 
be inferred by inspecting see Fig. [T] Surprisingly, we notice 
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Figure 9. Results for the general Alfven test, problem 4, at t = 
0.5 on 800 computational zones. The panels arc structured in a 
way similar to Fig. [5] Top panel: density, gas pressure and total 
pressure. Mid panel: x, y and z velocity components. Bottom 
panel: Lorentz factor 7 and transverse components of magnetic 
field. 



that error HLLD performs slightly better than HLLC. The 
numerical undershoots in density, in fact, are found to be 
~ 24% (HLLD) and ~ 32% (HLLC). The HLL solver is less 
prone to this pathology most likely because of the larger 
numerical diffusion, see the left panel close-up of Fig. [8] 

Errors (fo r B^) are computed using the exact solution 
available from ICiacomazzo fc RezzoUal (|2006l 'l which is free 
from the pathology just discussed. As shown in the bottom 
left panel of Fig. 1111 HLLD performs as the best numerical 
scheme yielding, at the largest resolution employed (3200 
zones), L-1 norm errors of ~ 18% to be compared to ~ 32% 
and ~ 46% of HLLC and HLL, respectively. 

The CPU times for the different solvers on this problem 
follow the proportion thii : ihiic : ^hiid = 1 : 1.1 : 1.4. 



4.1.5 Shock Tube 4 

The fourth sh ock tube test is taken f rom the "Generic 
Alfven" test in ICiacomazzo fc RezzoUal (|2006l ). The break- 
ing of the initial discontinuous states leads to the formation 
of seven waves. To the left of the contact discontinuity one 
has a fast rarefaction wave, followed by a rotational wave 
and a slow shock. Traveling to the right of the contact dis- 
continuity, one can find a slow shock, an Alfven wave and a 
fast shock. 

We plot, in Fig.|9l the results computed with the HLLD, 
HLLC and HLL Riemann solvers at t = 0.5, when the outer- 
most waves have almost left the outer boundaries. The cen- 
tral structure (0.4<2;<0.6) is characterized by slowly moving 
fronts with the rotational discontinuities propagating very 
close to the slow shocks. At the resolution employed (800 
zones), the rotational and slow modes appear to be visi- 




Figure 10. Magnification of the central region of Fig. |9] The 
left panel shows the density profile where the two slow shocks 
and the central contact wave are clearly visible. Central and right 
panels display the y components of velocity and magnetic field. 
Rotational modes can be most clearly distinguished only with the 
HLLD solver aX x 0.44 and x ^ 0.59. 







1/Ax 



Figure 11. L-1 norm errors (in 10^) for the four shock tube 
problems presented in the text as function of the grid resolution. 
The different combinations of lines and symbols refer to HLLD 
(solid, circles), HLLC (dashed, crosses) and HLL (dotted, plus 
signs). 



ble and distinct only with the HLLD solver, whereas they 
become barely discernible with the HLLC solver and com- 
pletely blend into a single wave using the HLL scheme. This 
is better shown in the enlargement of and profiles 
shown in Fig. 1101 rotational modes are captured at a:: ~ 0.44 
and X ~ 0.59 with the HLLD solver and gradually disappear 
when switching to the HLL scheme. 

At the contact wave HLLD and HLLC behave similarly 
but the sharper resolution attained at the left-going slow 
shock allows to better capture the constant density shell 
between the two fronts. 

Our scheme results in the smallest errors and numerical 
dissipation and exhibits a slightly faster convergence rate, 
see the plots in the bottom right panel of Fig. 1111 At low res- 
olution the errors obtained with HLL, HLLC and HLLD are 
in the ratio 1 : 0.75 : 0.45 while they become 1 : 0.6 : 0.27 as 
the mesh thickens. Correspondingly, the CPU running times 
for the three solvers at the resolution shown in Table |4] have 
found to scale as thii : thiic : ihiw = 1 : 1-4 : 1.8. This exam- 
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Figure 14. One dimensional cuts along the y (left) and z 
(right) axis showing the density profiles at different resolutions 
(128'^, 256'^ and 512'^) and with different solvers. Solid, dashed 
and dotted lines are used for the HLLD solver whereas plus and 
star symbols for HLL. 
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Figure 12. The 3D rotor test problem computed with HLLD 
(top panels) and HLL (bottom panels) at the resolution of 256^ . 
Panels on the left show the density map (at t = 0.4) in the xy 
plane at 2: = while panels to the right show the density in the 
xz plane at j; = 0. 
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Figure 13. Same as Fig. ll2l but showing the total pressure in the 
xy (left) and xz (right) panels for the HLLD solver. 



pie demonstrates the effectiveness and strength of adopting 
a more complete Riemann solver when describing the rich 
and complex features arising in relativistic magnetized flows. 



4.2 Multidimensional Tests 

We have implemented our 5 wave Riemann solver into the 
fram ework provided by the PLUTO code (|Mignone et al.l 
I2OO7I ). The constrained transport method is used to evolve 
the magnetic field. We use the third-order, total variation 
diminishing Runge Kutta scheme together with piecewise 
linear reconstruction. 



4.2.1 The 3D Rotor Problem 

We consider a three dimensional ver sion of the standard ro- 
tor problem (|Del Zanna et al1l2003l ). The initial condition 
consists of a sphere with radius ro = 0.1 centered at the 
origin of the domain taken to be the unit cube [—1/2, 1/2]^. 
The sphere is heavier (p — 10) than the surrounding {p — 1) 
and rapidly spins around the z axis with velocity compo- 
nents given by {v^ jV") = to {—y,x,0) where to — 9.95 
is the angular frequency of rotation. Pressure and magnetic 
field are constant everywhere, pg — 1, B — (1,0,0). 

Exploiting the point symmetry, we carried computa- 
tions until t = 0.4 at resolutions of 128^,256^ and 512^ 
using both the HLLD and HLL solvers. We point out that 
the HLLC of MB failed to pass this test, most likely because 
of the flux-singularity arising in 3D computations in the zero 
normal field limit. 

As the sphere starts rotating, torsional Alfven waves 
propagate outward carrying angular momentum to the sur- 
rounding medium. The spherical structure gets squeezed 
into a disk configuration in the equatorial plane {z = 0) 
where the two collapsing poles collide generating refiected 
shocks propagating vertically in the upper and lower half- 
planes. This is shown in the four panels in Fig. [12] show- 
ing the density map in the xy and xz planes obtained with 
HLLD and HLL and in Fig. [13] showing the total pressure. 
After the impact a hollow disk enclosed by a higher density 
shell at z = ±0.02 forms (top right panels in Fig I12p . In 
the xy plane, matter is pushed in a thin, octagonal-like shell 
enclosed by a tangential discontinuity and what seems to be 
a slow rarefaction. The whole configuration is embedded in 
a spherical fast rarefaction front expanding almost radially. 
Flow distortions triggered by the discretization on a Carte- 
sian grid are more pronounced with HLLD since we expect 
it to be more effective in the growth of small wavelength 
modes. 

In Fig. [14] we compare the density profiles on the y 
and 2 axis for different resolutions and schemes. From both 
profiles, one can see that the central region tends to be- 
come more depleted as the resolution increases. Inspecting 
the profiles in the y direction (left panel), we observe that 
HLL and HLLD tend to under- and over-estimate (respec- 
tively) the speed of the thin density shell when compared to 
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HLLD, t= 5. 


^^^^^^^^^^^^ 
^^^^^ 


^^^^^^^^^^^^^^ 


HLL, t= 5, 


HLL, t=15. 
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Figure 15. Color scale maps of B'^ + / Bz at different in- 
tegration times, 4 = 5, 15, 30. Panels on top (bottom) refer to 
computations accomplished with HLLD (HLL). Poloidal mag- 
netic field lines overlap. 



the reference solution computed with the HLLD solver at a 
resolution of 512^. The height of the shell peak is essentially 
the same for both solvers, regardless of the resolution. 

On the contrary, the right panel of Fig. [14] shows a sim- 
ilar comparison along the vertical z axis. At the same res- 
olution, HLL under-estimates the density peak located at 
z = 0.02 and almost twice the number of grid zones is needed 
to match the results obtained with the HLLD solver. The 
location of the front is approximately the same regardless of 
the solver. 

In terms of computational cost, integration carried with 
the HLLD solver took approximately ~ 1.6 that of HLL. 
This has to be compared with the CPU time required by 
HLL to reach a comparable level of accuracy which, dou- 
bling the resolution, would result in a computation ~ 2* 
as long. In this respect, three dimensional problems like the 
one considered here may prove specially helpful in establish- 
ing the trade off between numerical efficiency and accuracy 
which, among other things, demand choosing between ac- 
curate (but expensive) solvers versus more diffusive (cheap) 
schemes. 



4-2.2 Kelvin- Helmholtz Unstable Flows 

The setup, taken from lBucciantini fc Del Zann al (|2006h . con- 
sists of a 2D planar Cartesian domain, x £ [0, 1], y G [^1,1] 
with a shear velocity profile given by 



0.01 



10" 
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HLLD(H 
HLLD(M) 
HLLD(L) 




Figure 16. Top: growth rate (as function of time) for the Kelvin- 
Helmholtz test problem computed as Av^ = (t^max — ■ )/2 at 
low (L), medium (M) and high (H) resolutions. Solid, dashed and 
dotted lines show results pertaining to HLLD, whereas symbols 
to HLL. Bottom: small scale power as a function of time for the 
Kelvin-Helmholtz application test. Integrated power is given by 
Ps = 1/2 ° 1^ \ l^e^' y)\^'^y'^^ where V{k, y) is the complex, 
discrete Fourier transform of v^^x, y) taken along the x direction. 
Here ks is the Nyquist critical frequency. 



Density and pressure are set constant everywhere and ini- 
tialized to p = 1, pg = 20, while magnetic field components 
are given in terms of the poloidal and toroidal magnetization 
parameters o-poi and ator as 



(B",B'> 



(61) 



where we use a^oi ~ 0.01, ator = 1. The shear layer is per- 
turbed by a nonzero component of the velocity. 



— sin (27ra;) exp 

400 ^ ' ^ 



y_ 

>2 



(62) 



v"^ = --tanh(lOOy) 



(60) 



with /3 — 1/10, while we set — 0. Computations are 
carried at low (L, 90 x 180 zones), medium (M, 180 x 360 
zones) and high (H, 360 x 720 zones) resolution. 

For t<5 the perturbation follows a linear growth phase 
leading to the formation of a multiple vortex structure. 
In the high resolution (H) case, shown in Fig 1151 we ob- 
serve the formation of a central vortex and two neigh- 
bor, more stretched ones. These elongated vortices are not 
seen in the computation of iBucciantini fc Del Zannal (|2006l ) 
who employed the HLL solver at our medium resolution. 
As expected, small scale patterns are best spotted with 
the HLLD solver, while tend to be more diffused using 
the two-wave HLL scheme. The growth rate (computed as 
Av^ = (t^max — *'min)/2i See top panel in Fig. \T^ . is closely 
related to the poloidal field amplification which in turn pro- 
ceeds faster for smaller numerical resistivity (see the small 
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sub-plot in the same panel) and thus for finer grids. Still, 
computations carried with the HLLD solver at low (L), 
medium (M) and high (H) resolutions reveal surprisingly 
similar growth rates and reach the saturation phase at es- 
sentially the same time (t 3.5). On the contrary, the sat- 
uration phase and the growth rate during the linear phase 
change with resolution when the HLL scheme is employed. 

Field amplification is prevented by reconnection events 
during which the field wounds up and becomes twisted by 
turbulent dynamics. Throughout the saturation phase (mid 
and right panel in Fig llSp the mixing layer enlarges and the 
field lines thicken into filamentary structures. Small scale 
structure can be quantified by considering the power residing 
at large wave numbers in the discrete Fourier transform of 
any fiow quantity (we consider the y component of velocity). 
This is shown in the bottom panel of Fig [16] where we plot 
the integrated power between ks/2 and ks as function of 
time {ks is the Nyquist critical frequency). Indeed, during 
the statistically steady flow regime {t>20), the two solvers 
exhibits small scale power that differ by more than one order 
of magnitude, with HLLD being in excess of lO"'^ (at aU 
resolutions) whereas HLL below 10~^. 

In terms of CPU time, computations carried out with 
HLLD (at medium resolution) were ~ 1.9 slower than HLL. 



4-2.3 Axisymmetric Jet Propagation 

As a final example, we consider the propagation of a rela- 
tivistic magnetized jet. For illustrative purposes, we restrict 
our attention to axisymmetric coordinates with r G [0, 20] 
and 2 G [0, 50]. The jet initially fills the region r,z ^1 with 
density pj — 1 and longitudinal {z) velocity specified by 

-yj = 10 (■u" = = 0). 

The magnetic field topology is described by a constant 
poloidal term, B" , threading both the jet and the ambient 
medium and by a toroidal component B'''(r) = 'yjb^(r) with 



0.0 0.2 0,3 0.5 0,7 0.8 1.0 




for 
for 



r < a , 

a <r <1, 



(63) 



where a = 0.5 is the magnetization radius and bm is a con- 
stant and vanishes outside the nozzle. The thermal pressure 
distribution inside the jet is set by the radial momentum 
balance, rdrPg = —b^dr{rb^) yielding 



Pair) =pj+h„ 



1 — min — , 1 



(64) 



where pj is the jet/ambient pressure at r = 1 and is re- 
cove red from the definit ion of the Mach number, M = 
'^jVpj/(rpj) + l/(r-l), with M = 6 and r = 5/3, al- 
though we evolve the equa tions using the approximate d 
Synge gas equation of state o f lMignone fc McKinnevI l|2007h . 

The relative contribution of the two components is 
quantified by the two average magnetization parameters 
a. = Bl/{2 (p,)) = (hi) /(2 (p,)) yielding 



6m = 



a2(2cr^ - 1 +41oga) 



= V'^z {bl,a^ + 2pj) , (65) 



where for any quantity q{r), (q) gives the average over the 
jet beam r G [0,1]. We choose cr^ = 0.3, = 0.7, thus 
corresponding to a jet close to equipartition. 




-40 - 



Figure 17. Left: composite color map image of the jet at t = 270 
at the resolution of 40 points per beam radius. In clockwise direc- 
tion, starting from the top right quadrant: density logarithm, gas 
pressure logarithm, thermal to total pressure ratio and <j> compo- 
nent of magnetic field. The color scale has been normalized such 
that the maximum and minimum values reported in each subplots 
correspond to 1 and 0. 



The external environment is initially static {ve = 0), 
heavier with density pe = 10^ and threaded only by the 
constant longitudinal field B^. Pressure is set everywhere to 
the constant value pj . 

We carry out computations at the resolutions of 10, 20 
and 40 zones per beam radius (r = 1) and follow the evo- 
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Figure 18. Enlargement of the turbulent flow region [2, 10] X 
[10, 18] at t = 300 showing the poloidal magnetic field structure 
(in log scale) for the high and medium resolution runs (40 and 20 
points per beam radius. 




50 100 150 200 250 300 
time 

Figure 19. Volume average of VB^/Bp as a function of time. 
Here Bp is the poloidal magnetic field. Solid, dashed and dot- 
ted lines refers to computations carried out with HLLD, whereas 
symbols give the corresponding results obtained with HLL. 



lution until t — 300. The snapshot in Fig. [17] shows the 
solution computed ai t — 300 at the highest resolution. 

The morphological structure is appreciably affected by 
the magnetic field topology and by the ratio of the mag- 
netic energy density to the rest mass, b^/p « 0.026. The 
presence of a moderately larger poloidal component and a 
small Poynting flux favor the formation of a hammer-like 



structure rather than a nose cone (see lLeismann et al.ll2005l : 
iMignone et al.l [20051 ). At the termination point, located at 
2 ~ 40.5, the beam strongly decelerates and expands radi- 
ally promoting vortex emission at the head of the jet. 

Close to the axis, the flow remains well coUimated 
and undergoes a series of deceleration/acceleration events 
through a series of conical shocks, visible at z ~ 
4.5, 19, 24, 28, 32. Behind these recoUimation shocks, the 
beam strongly decelerates and magnetic tension promotes 
sideways defiection of shocked material into the cocoon. 

The ratio Pg/p (bottom left quadrant in Fig [17]) clearly 
marks the Kelvin-Helmholtz unstable slip surface separat- 
ing the backflowing, magnetized beam material from the 
high temperature (thermally dominated) shocked ambient 
medium. In the magnetically dominated region turbulence 
dissipate magnetic energy down to smaller scales and mixing 
occurs. The structure of the contact discontinuity observed 
in the figures does not show suppression of KH instabil- 
ity. This is likely due to the larger g rowth of the toroida l 
field component over the poloidal one (jKeppens et al.ir2008l ). 
However we also think that the small density ratio (10~^) 
may favor the growth of instability and moment um transfer 
throu gh entrainment of the external medium (jRossi et al.l 
l2008l ). 

For the sake of comparison, we also plot (Fig I18|l the 
magnitude of the poloidal magnetic field in the region r G 
[2, 10], 2 G [10, 18] where turbulent patterns have developed. 
At the resolution of 40 points per beam radius, HLLD dis- 
closes the finest level of small scale structure, whereas HLL 
needs approximately twice the resolution to produce similar 
patterns. This behaviour is quantitatively expressed, in Fig 
1191 by averaging the gradient log(_B^ + _B?) over the volume. 
Roughly speaking, HLL requires a resolution ~ 1.5 that of 
HLLD to produce pattern with similar results. 



5 CONCLUSIONS 

A five-wave HLLD Riemann solver for the equations of rel- 
ativistic magnetohydrodynamics has been presented. The 
solver approximates the structure of the Riemann fan by 
including fast shocks, rotational modes and the contact dis- 
continuity in the solution. The gain in accuracy comes at the 
computational cost of solving a nonlinear scalar equation in 
the total pressure. As such, it better approximates Alfven 
waves and we also found it to better capture slow shocks 
and compound waves. The performance of the new solver 
has been tested against selected one dimensional problems, 
showing better accuracy and convergence properties than 
previously known schemes such as HLL or HLLC. 

Applications to multi-dimensional problems have been 
presented as well. The selected tests disclose better resolu- 
tion of small scale structures together with reduced depen- 
dency on grid resolution. We argue that three dimensional 
computations may actually benefit from the application of 
the proposed solver which, albeit more computationally in- 
tensive than HLL, still allows to recover comparable accu- 
racy and resolution with a reduced number of grid zones. 
Indeed, since a relative change S in the mesh spacing results 
in a factor 5* in terms of CPU time, this may largely favour 
a more sophisticated solver over an approximated one. This 
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issue, however, need to receive more attention in forthcom- 
ing studies. 
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APPENDIX A: PROPAGATION OF 
ROTATIONAL DISCONTINUITIES 

Left and right states across a rotational discontinuity can be 
found using the results outlined in i]3.2l More specifically, we 
construct a family of solutions parameterized by the speed of 
the discontinuity and one component of the tangential 
field on the right of the discontinuity. Our procedure ca n 
be shown to be be equivalent to that of I Komissarov! (|l997h . 
Specifically, one starts by assigning p,pg,v, on the left 
side of the front (B* = (0, , B^)) together with the speed 
of the front, . Note that cannot be freely assigned but 
must be determined consistently from Eq. (|46|) . Expressing 
K'' {k / x) in terms of v'' , B*" and and substituting back 
in the x— component of H46p . one finds that there are two 
possible values of B^ satisfying the quadratic equation 

a(B")' + feB" -f c = , (Al) 



where the coefficients of the parabola are 



2x ( 



V ■ 
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■ w„ + 



B* B' 
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(A2) 



(A3) 



with rj = 1- {v^f ~ {v^f, X = v^B^ + w^B^ and 7 being 
the Lorentz factor. The transverse components of K are 
computed as 



(A4) 



On the right side of the front, one has that p, pg, w, B^ 
and K are the same, see jj3.2l Since th e transverse field is 
elliptically polarized (jKomissarovl 1 199*71 ) . there are in prin- 
ciple infinite many solutions and one has the freedom to 
specify, for instance, one component of the field (say B^). 
The velocity vr and the z component of the field can be 
determined in the following way. First, use Equation (|47|l to 
express (k = x, y, z) as function of BJ for given Bfi and 
B]^. Using the jump condition for the density together with 
the fact that p is invariant, we solve the nonlinear equation 

/5L7i (K'' - Vl) = PR-JR (A"" - Vr) , (A5) 

whose roots gives the desired value of B|j. 



